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20.  ABSTRACT  (Continued) 

i 

I 

central  weighted  first-step  Lax-Wendroff  time  integration  with  parabolic  spatial  discretization  either 
in  its  full  or  condensed  [M]  matrix  form.  Both  the  standard  and  central  weighted  first-step 
Godunov  time  integrations  are  found  to  be  numerically  diffusive.  This  diffusivity  tends  to 
override  whatever  spatial  discretization  is  used.  However,  the  positivity  property  possessed  by  the 
Godunov  schemes  can  be  valuable  for  many  applications. 
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A  COMPARISON  OF  TWO-STEP  TIME 
INTEGRATION  SCHEMES  FOR  THE  FINITE  ELEMENT 
ADVECTION  EQUATION 


1.  INTRODUCTION 

Nonlinear  hyperbolic  equations  of  the  form 

+div(Fp)=  ?(K,p)  (1.1) 

ot 

where  p  is  a  vector  of  conserved  quantities,  Kis  the  velocity  field,  and  Fis  a  given  functional,  describe 
flow  behavior  in  subjects  ranging  from  hydraulics  to  gas  dynamics.  Finite  difference  techniques  that 
had  been  applied  to  the  solution  of  Equation  (1.1)  were  examined  in  detail  by  Roache  [1]  in  his  classic 
text  Computational  Fluid  Dynamics  published  in  1972.  More  recent  finite  difference  attacks  on  Equation 
(1.1)  have  been  reviewed  by  Sod  [2]  and  by  Book,  et  al.  [3].  As  is  evident  from  these  references,  finite 
difference  solutions  of  Equation  (1.1)  have  reached  a  high  level  of  sophistication  and  accuracy.  The 
most  significant  criticism  that  can  be  leveled  against  these  techniques  is  the  difficulty  encountered  in 
treating  complex  flow  boundaries.  This  is  especially  true  for  the  most  accurate  methods  that  rely  on  the 
use  of  a  "staggered"  mesh. 

A  natural  way  of  handling  complex  flow  boundaries  is  through  finite  element  discretization  of  the 
spatial  derivatives  in  Equation  (1.1).  In  these  finite  element  methods  (for  an  excellent  introduction  see 
the  text  by  Baker  [4])  the  space  of  concern  is  divided  into  a  large  number  of  sub-spaces  over  each  of 
which  the  dependent  flow  variables  are  approximated  by  shape  functions.  Complex  flow  boundaries  fall 
naturally  out  of  the  discretization  as  constraints  on  the  shape  function  coefficients  of  boundary  ele¬ 
ments. 
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Before  blindly  applying  finite  element  methods  to  the  solution  of  Equation  (1.1),  however,  it  is 
necessary  to  evaluate  their  ability  to  replicate  some  simple  flow  behaviors  predicted  by  Equation  (1.1). 
To  this  end,  we  examine  in  this  report  numerical  solutions  of  the  simple,  constant  velocity,  advection 
equation 

+  V  I2-  =  0  (1.2) 

at  ax 

that  becomes,  upon  finite  element  spatial  discretization, 

[M]  («]  +  V  [AT]  {/*)  =  0.  (1.3) 

Here,  [Ml  and  [K]  are  matrices  that  depend  on  the  specific  discretization,  (/?[  is  a  column  matrix 

representing  the  node  point  values  of  the  "density"  p,  and  the  dot  denotes  ordinary  differentiation  with 

respect  to  time.  For  practical  computational  purposes,  two  requirements  can  be  imposed  on  the  time 

integration  technique  used  to  advance  Equation  (1.3): 

1.  It  should  require  knowledge  of  (ftj  only  at  the  node  points.  This  requirement  is  imposed  to 
preserve  the  strongest  attributes  of  the  finite  element  method,  namely,  ease  of  establishing 
grids  for  complex  geometric  boundaries  and  ease  of  handling  boundary  conditions. 

2  It  should  be  explicit  to  ensure  fast,  efficient  calculations.  This  requirement  really  stems  from 
Equation  (1.1)  where  an  implicit  integration  scheme  would  mandate  the  solution  of  large, 
nonlinear  sets  of  algebraic  equations  at  each  time  step. 

The  problems  of  concern  are  a)  what  combination  of  spatial  discretization  and  time  integration 
"best"  models  the  true  advection  solution  and  b)  whether  "best"  is  good  enough.  Within  bounds,  the 
latter  concern  is,  of  course,  fairly  subjective. 

In  the  subsequent  sections  of  this  report,  we  undertake  numerical  experiments  using  both  linear 
and  parabolic  spatial  discretizations  of  Equation  (1.2)  combined  with  two  basic,  two-step  time  integra¬ 
tion  schemes  for  Equation  (1.3).  One  of  these  schemes  is  first  order  accurate  (in  a  Taylor  series  sense) 
while  the  other  is  second  order  accurate.  Following  the  convention  of  Sod  (2).  we  refer  to  the  first 
scheme  as  a  Godunov  type  scheme  and  the  second  scheme  as  a  Lax-Wendroff  type  scheme  We  also 
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consider  the  effects  of  [AT]  matrix  condensation  [4]  on  the  numerical  solution.  This  popularly  used 
technique  obviates  the  need  for  performing  an  inversion  of  the  large  [AT]  matrix.  In  addition,  two  first 
step  weightings  of  the  basic  time  integration  schemes  are  studied.  Such  weightings  are  used  commonly 
in  finite  difference  methods  [1]  to  average  the  influence  of  adjacent  nodes.  One  weighting  gives  zero 
weight  to  the  node  under  consideration  and  is  termed  the  "smoothed"  weighting.  The  second  weight¬ 
ing,  based  on  symmetry  considerations,  is  termed  the  "modified"  weighting 

All  told,  twenty  four  different  combinations  of  spatial  discretization  and  time  integration  are 
applied  to  the  advection  of  both  a  smooth  cosine  hill  density  distribution  and  a  discontinuous  square 
hill  density  distribution.  For  the  most  part,  those  combinations  that  employ  a  parabolic  spatial  discreti¬ 
zation  yield  more  accurate  solutions  than  those  that  employ  a  linear  spatial  discretization.  This  result  is 
similar  to  some  recent  findings  of  Leonard  [5]  that  indicate  even-ordered  finite  differencing  to  be  more 
accurate  than  odd-ordered  finite  differencing.  Overall,  the  "best”  representation  of  the  true  advection 
solution  is  obtained  from  the  combination  of  parabolic  spatial  discretization  and  modified  Lax-Wendroff 
lime  integration  This  time  integration  combined  with  the  condensed  parabolic  I  A/]  matrix  formulation 
also  gives  quite  acceptable  numerical  solutions.  The  modified  Godunov  time  integration  combined  with 
either  the  full  or  condensed  parabolic  spatial  discretization,  while  highly  diffusive,  does  possess  some 
positivity  properties  that  can  be  useful  in  certain  applications. 

We  conclude  the  report  by  examining  time  step  (Courant  number)  stability  questions  for  the 
more  promising  combinations  of  spatial  discretization  and  time  integration. 

2.  PROBLEM  DEFINITION 

Let  the  density  of  a  certain  substance  be  denoted  by  />(x,r)  This  substance  is  being  advected 
downstream  ( v-direction)  at  a  constant  velocity  V.  Requiring  that  the  substance  be  conserved  leads  to 
the  advection  equation: 

&-+V&--0.  (2.1) 

dt  di 
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Normally  this  equation  is  applied  to  a  finite  (or  infinite)  length  in  the  x-direction,  thus  permitting 
mass  to  flow  out  at  the  boundaries.  This  would  confuse  the  evaluation  of  the  integration  schemes,  in 
that  the  resulting  errors  could  either  be  in  the  scheme  itself  or  the  system  could  be  losing  mass.  In 
order  to  avoid  the  possibility  of  gaining  or  losing  mass,  the  x  domain  of  Equation  (2.1)  is  taken  as  a 
continuous  loop. 


Consider  now  that  the  x-direction  is  a  continuous  loop,  hereafter  called  a  racetrack.  The  domain 
of  x  is  0  <  x  <  L.  By  this  we  mean  that  the  point  x  =  L  is  the  same  as  the  point  x  =  0  as  shown  in 
Figure  2.1.  For  our  subsequent  numerical  studies,  we  take  L  =  48  (arbitrary  units  of  distance). 


Two  initial  conditions  are  considered:  (a)  a  cosine  hill  between  x,  and  x2, 


p  (x,  0)  =  F(x) 


1  x  ^  X[  and  x  ^  x2 

3  i  27r  (x  -  X)) 

—  —  —  COS - X j  <  X  <  x2 

2  2  x2  —  x  | 


and  (b)  a  square  hill  also  between  xx  and  x2, 

1  x  <  X!  and  x  ^  x2 
p  (x,  0)  =  Fix)  =  j2  Xi<Jf<X2 

These  initial  conditions  are  shown  in  Figure  2.2  for  X|  =  8  and  x2  =  18. 


(2.2) 


(2.3) 


The  exact  solution  to  Equation  (2.1)  is  expressed  as 

p(x,f)  =  F(x  -  Vt)  (2.4) 

where  F  is  the  initial  density  distribution  given  by  Equation  (2.2)  or  (2.3).  Note  that  the  values  of  x  in 

Equation  (2.4)  must  reflect  the  racetrack  domain  and  the  racetrack  boundary  condition 

p(0.r)  =  p(Z..r).  (2.5) 

Although  the  exact  solution  is  near  to  being  trival,  the  numerical  solution  is  extremely  challeng¬ 
ing  as  witnessed  by  the  literature  (see  References). 


3.  FINITE  FI  TMENT  FORMULATION 


The  Finite  I  lenient  Method  is  a  technique  which  divides  the  space  domain  into  subdomains. 


called  elements  Over  each  element  the  density  is  approximated  by 
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Pe(x,t)=  £  Rf(t)  JV,(x)  (3.1) 

f“  I 

where  pe(x,t)  is  the  approximation  of  p(x,t)  over  the  element  e ,  Rf(t)  are  a  set  of  discrete  values  at 
nodes  (or  grid  points)  on  the  boundaries  or  within  the  element,  and  N,(x)  (which  are  taken  to  be  the 
same  form  for  each  element)  are  a  set  of  interpolation  polynomials,  also  called  shape  functions,  of 
order  K  —  1  where  A  is  the  number  of  nodes  per  element.  The  following  subsections  will  define  these 
symbols  for  each  shape  function  assumed;  i.e. ,  linear  and  parabolic. 

In  general  the  approximate  solution  will  be  in  error  and  will  not  satisfy  Equation  (2.1).  Thus  we 
can  write, 

+  V  =  E‘lx.1).  (3.2) 

Ot  ox 

For  any  assumed  p‘  we  seek  a  solution  for  R which  will  minimize  E‘(x,t).  There  are  several  tech¬ 
niques  to  do  this;  however.  Galerkin’s  method  is  used  here.  The  method  requires  that  the  error 
E'(x,t )  be  orthogonal  with  iV,(x)  for  each  element;  therefore, 

f '  E'(x.i)  )V,(x)  dx  =  0  /  =  1 . K  (3.3) 

•Jii  ' 

where  L''  is  the  element  length.  Substitution  of  Equation  (3.1)  into  (3  3)  gives  the  classic  Finite  Ele¬ 


ment  discretization  equation 

l  AM  {/?'}  +  V  (A'  l  | AT)  =  (0)  G  4) 

where 

I  AM  —  <  =  ,V(  (a  )  iV,(x)  dx  i.j  =  1 . K  (3.5  > 

[A1  —  k‘„  =  f  N‘  (x)  .V,  <x>  dx  i.j  =  I . A  (3.6) 


(A'l  =  (A|  R'2  . R'kV 

(A'l  =  (A'i  R)  . A*)' 

The  prime  denotes  spatial  differentiation  and  the  dot  denotes  ordinary  time  differentiation 
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When  equation  (3.4)  is  written  for  each  element  and  connectivity  between  elements  is  considered 
then  the  global  system  of  equations  to  be  solved  is 

[A#]  l«l  +  V\K\  (/?)  =  (0)  (3.7a) 

where 

r 

{/?!  =  R{R2  ....  (3.7b) 

Here  [M 1  and  [ A" ]  are  the  assembled  mass  and  advection  matrices  and  N  is  the  total  number  of  nodal 
values  (degrees  of  freedom)  in  the  system. 

3.1  Linear  Elements 


Linear  finite  elements  are  developed  by  using  a  linear  shape  function  with  the  coefficients 
expressed  in  terms  of  nodal  values.  Thus  for  the  linear  element  we  have  from  Equation  (3.1)  with 
A  =  2: 

fi'  (x.  t)  —  R\  ,V,  +  R':  Ni 

where 

.V,  =  1  -  (x/L)  (3.8) 

,V,=  (x/Lr) 

and  where  R\  is  the  density  at  ,v  =  0  and  R':  is  the  density  at  v  =  L‘\  (Again,  Ll  is  the  element 
length).  Integrating  Equation  (3.5)  with  the  shape  function  of  Equation  (3.8)  defines  the  mass  matrix 
for  i.i  =  1  to  2 
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The  assembled  mass  and  advection  matrices  appearing  in  Equation  (3.7)  are  given,  for  the  linear  ele¬ 
ment  discretization,  in  Appendix  A. 

3.2  Parabolic  Elements 

Parabolic  finite  elements  are  developed  by  using  parabolic  shape  functions  with  coefficients 
expressed  in  terms  of  nodal  values.  For  this  element  a  node  at  the  mid-point  is  established  in  addition 
to  the  nodes  at  the  ends.  Equation  (3.1)  is  written  for  K  =  3: 

p'lx.f)  =  R\  N |  +  R‘,  N j  +  R\  /V, 

where 

A,  =  1  -  3 (x/L1)  +  2 (x/L1)2 

N2  =  4(x/L‘)  -  4(x/L')2  (3.12) 

N}~-  (x/L1)  +  2 (x/Lc)2 

and  where  R\,R'2  and  R\  are  the  densities  at  a  =  0.  a  =  U/ 2,  and  x  —  L'\  respectively.  Integrating 
Equation  (3.5)  with  the  above  shape  functions  defines  the  mass  matrix  for  /,/  =  1,2.3: 

/.  4  2  -1 

[M'1  =  —  2  16  2  (3.13) 

30  l_  i  2  4 

Likewise,  integrating  Equation  (3.6)  after  differentiating  Equations  (3.12)  defines  the  advection  matrix 
for  ij  =  1,2,3: 

-.1-3  4  -1 

1A'1]  =  —  -4  0  4  .  (3.14) 

6  l  1  -4  3 

Now  Equation  (3.4)  can  be  written  for  the  parabolic  element: 


2  -l]  R‘'  v  -3  4  -l]  R'\  0 

16  2  A3+  —  -4  0  4  R\  =  O'.  (3  IS) 

2  4  R\  6  ,  _4  3  R,  o 

The  assembled  mass  and  advection  matrices  appearing  in  liquation  (3.7)  are  given,  for  the  parabolic 
element  discretization,  in  Appendix  B. 

3.3  Condensed  I  3/ 1  Matrix  Formulation 


Equation  <3  7  I  can  he  rew  ritten  lormalU  as 


i 

\ 


R 


‘  [Ml  1 1 A  1 1  A  | 


(316) 
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where  [A/1-1  denotes  the  inverse  of  (A/1.  Generally,  [A/1  is  a  large,  sparse  matrix  and  the  calculation 
of  [A/] 1  (which  is  a  full  matrix)  can  be  very  time  consuming.  The  process  of  condensing  the  mass 
matrix,  or  lumping  the  mass  at  the  diagonal  terms,  is  a  popular  procedure  [4]  for  avoiding  this  time 
consuming  inversion. 


Hence,  the  condensed  mass  matrix  [A/,]  is  defined  by 


[A/, 


w,  „  =  L  >w„  i  -  1 . N 

mul  =  0.  i  j*  j  i.j  —  1 .  N 


(3.17) 


where  the  mti  are  the  entries  in  [A/]  and  N  is  the  rank  of  [A/1.  The  matrix  [Mc\  now  replaces  [A/]  in 
Equation  (3.7)  and  Equation  (3.16)  is  replaced  by 


\R)  =  -  V[Mc]-'[K)[R).  (3.18) 

We  note  that  this  condensation  is  mass  preserving. 


4.  TIME  INTEGRATION  SCHEMES 


There  are  many  methods  available  for  advancing  Equation  (3.16)  [or  (3.18)1  in  time.  For  reasons 
explained  in  the  Introduction,  the  focus  of  this  report  deals  with  the  two-step  methods  of  Lax-Wendroff 
and  Godunov  (see  Reference  [2]). 


4.1  Standard  Two-Step 


Consider  that  the  value  of  the  nodal  densities  at  time  step  n  are  known  {/?"}.  The  solution  at  the 
n  +  1  time  step  is  sought.  The  time  increment  is  Af.  Then  the  standard  two-step  method  can  be  writ¬ 
ten  as: 


1st  step 

{R"H'} 

=  (/?")  +a  \l{R") 

(4.1a) 

2nd  step 

=  1/?")  +  Af|K',+,,l 

(4.1b) 

where  «  is  a  fraction  of  the  time  step.  The  time  derivative  \R")  is  determined  from  Equation  (3.16) 


[or  (3.18)1  using  \R’’}  and  subsequently  (/?"*")  is  determined  from  Equat’on  (3.16)  [or  (3.18)1  using 


I/?""’)  from  Equation  (4  la). 
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The  two  well-known  methods  embodied  in  equations  (4.1)  are  the  Lax-Wendroff  method  when 
a  ~  1/-  and  the  Godunov  method  when  a  =  1.  The  Lax-Wendroff  method  first  estimates  the  time 
derivative  at  the  mid-point  of  the  time  step  then  uses  the  value  to  step  ahead  to  n  +  1.  The  Godunov 
method  first  estimates  the  time  derivative  at  a  full  step  and  uses  this  estimate  to  step  ahead  to  n  +  1. 

Equations  (4.1)  are  referred  to  in  this  report  as  the  standard  two-step  methods.  Because  the  stan¬ 
dard  two-step  methods  update  a  nodal  value  only  as  a  function  of  its  previous  value,  it  seems  plausible 
that  introducing  a  weighting  matrix  to  average  the  influence  of  adjacent  nodes  might  improve  numerical 
accuracy.  Such  averaging  is  used  frequently  inite  difference  methods  [1).  This  weighting  matrix 
should  have  a  general  form  so  as  to  be  problem  independent.  Two  such  weighting  matrices  are 
employed  herein. 

4.2  Smoothed  Two-Step 

The  first  weighting  matrix  is  termed  the  "smoothed"  weighting  matrix.  The  terms  in  the 
smoothed  weighting  matrix,  [  Ws\,  are  determined  from  the  entries  in  the  global  mass  matrix  [A/]: 

ww  =  0  t'=l . N 

m,, 

ws"  =  i  *  j  i,  j  —  1 . N 

C  =  £  m,r  (4.2) 

7-1 
/  ^  I 

The  quantity  C,  is  the  sum  of  all  terms  in  row  'i’  of  [M]  except  the  diagonal  term.  Note  that  the  sum 
across  row  ‘i’  of  [  W/v]  is  unity  so  that  [  W'',]  is  mass  conserving. 

The  weighting  matrix  is  used  onlv  in  the  first  step.  Thus  the  smoothed  two-step  method  is  written 
as: 

1st  step  (K”*")  =  [«']{/?'’}  +a  Ar{/T) 
and 

2nd  step  {«n  +  l}  =  {/?"(  +  At (4.3) 
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The  smoothed  Lax-Wendroff  method  is  obtained  by  letting  «  =  1/2  in  Equation  (4  3),  and  the 
smoothed  Godunov  method  is  obtained  by  letting  a  =  1  in  Equation  (4.3). 

4.3  Modified  Two  Step 

The  second  weighting  matrix  is  termed  the  "modified"  weighting  matrix.  The  terms  in  the 
modified  weighting  matrix  [  IT',,,],  are  determined  from  [A/]:  the  global  mass  matrix: 

(^,,1  -  =  ~~  U=  1 . A/ 

I 

V 

C,  =  £  mT  (4 .4) 

/-i 

Now  C  includes  the  diagonal  term  m„  and  the  sum  across  row  ‘i’  of  [  is  again  unity. 

The  weighted  matrix  is  again  used  only  in  the  first  step.  Thus  the  modified  two-step  method  is 
written  as. 

1st  step  I R . 1  =  liE'JlK"}  +  «  A/lK") 

and 

2nd  step  {R'"')=  {/?"}  +  (4.5) 

The  modified  Lax-Wendroff  method  is  obtained  by  letting  a  =  1/2  in  Equation  (4.5):  and  the 
modified  Godunov  method  is  obtained  by  letting  a  =  1  in  Equation  (4.5). 

5.  NUMERICAL  RESULT 

Numerical  computations  were  performed  to  evaluate  the  ability  of  the  various  solution  schemes  to 
replicate  the  true  advection  solution.  There  are  several  measures  by  which  to  ascertain  the  adequacy  of 
a  particular  scheme,  one  being  the  eyeball.  A  second  measure  used  herein  is  the  average  absolute 
error;  hereafter  called  simply  the  Error  given  by 

Error  =  -^  Ip  -  pi  (5.1) 

where  />  is  the  exact  solution,  p  is  the  approximate  solution,  and  N  is  the  number  of  nodes.  This  quan¬ 
titative  measure  is  used  to  compare  the  schemes.  The  scheme  with  the  smallest  Error  is  the  "best" 
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numerical  model.  At  which  value  this  error  becomes  intolerable  is  still  subjective.  However  by  observ¬ 
ing  the  density  distributions  at  a  given  time  and  the  associated  hr  tot.  an  acceptable  range  can  be  esta¬ 
blished. 

As  mentioned  previously,  the  length  of  the  racetrack  was  set  at  L  =  48  (arbitrary  units  of  dis¬ 
tance)  and  the  initial  location  of  the  cosine  hill  and  square  hill  density  distributions  [Equations  (2.21 
and  (2.3)]  was  between  =  8  and  v:  =  18.  The  advection  velocity  was  set  at  F  =  1  unit  of 
distance/s.  Hence,  the  density  distribution  is  advected  once  around  the  track  every  48  s. 

In  Section  5.1,  the  results  obtained  using  the  linear  finite  element  spatial  discretization  are 
presented.  Section  5  2  contains  the  results  obtained  using  the  parabolic  finite  element  spatial  discretiza¬ 
tion.  These  results  are  based  on  a  time  period  of  96  s  (two  cycles  around  the  track)  and  are  compared 
in  Section  5.3.  The  effect  of  time  step  size  on  the  better  solution  schemes  is  studied  in  Section  5.4. 
Finally,  Section  5.5  presents  results  for  long  time  periods.  The  density  distribution  is  allowed  to  travel 
around  the  track  for  480  s  (ten  cycles). 

5.1  Linear  Elements 

The  development  of  the  linear  finite  element  (linear  FEM)  spatial  discretization  is  given  in  Sec¬ 
tion  3.1  and  Appendix  A  We  take  the  length  of  each  element  as  L' -  1  unit  of  distance.  Thus,  there 
are  48  elements  on  the  track.  Node  1  corresponds  to  x  —  0  while  node  ,V  +  1  =  49  corresponds  to 
x  =  L  Since  the  track  is  closed,  the  number  of  degrees  of  freedom  N  =  48. 

Based  on  some  trial  calculations,  a  time  step  Ar  =  0.2  s  was  selected  to  study  the  time  integration 
schemes.  The  Courant  number,  equal  to  FA i/ L'\  is  0.2  for  the  linear  element. 

Figures  5.1.1  to  5.1.6  present  results  of  the  advected  cosine  hill  density  distribution  for  the  three 
Lax-Wendroff  methods  (standard,  smoothed,  and  modified)  and  the  three  Godunov  methods  using  the 
linear  FEM.  Figures  5.1.7  to  5.1.12  contain  similar  results  using  the  condensed  linear  mass  matrix  for¬ 
mulation  (linear  CFM).  It  is  worthwhile  to  point  out  that  the  linear  CFM  spatial  discretization  is 
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equivalent  to  standard  central  differencing.  The  time  periods  in  each  figure  are  24  s,  48  s,  72  s,  and  96 
s  or  1/2,  1,  1-1/2,  and  2  cycles  around  the  track.  The  exact  solution  is  superimposed  on  the  graphs  of 
the  advected  hill.  Using  a  time  step  of  0.2  s  means  at  96  s  the  solution  has  gone  through  480  time 
steps. 


The  first  item  to  notice  from  Figure  5.1.1  is  that  the  standard  Lax-WendrofT  method  does  an 
excellent  job  of  advecting  the  cosine  hill.  From  Figure  5.1.13  (a)  the  Error  at  96  s  is  1.4%.  Looking  at 
the  other  linear  FEM  figures,  the  remaining  time  integrators  tend  to  diffuse  the  hill.  This  is  numerical 
diffusion  due  to  the  nature  of  the  two-step  method.  The  standard  Godunov  scheme  is  the  only  one 
which  shows  no  oscillations.  Both  the  standard  Lax-Wendroff  scheme  and  standard  Godunov  scheme 
"track"  the  hill  (do  not  lead  or  lag),  but  the  latter  has  a  50%  reduction  of  the  peak  value  at  96  s.  Fig¬ 
ures  5.1.13  (a)  and  (b)  show  that  the  smoothed  integrator  is  the  poorest,  the  standard  is  the  best,  and 
the  modified  is  in  between.  Using  the  Error  and  the  density  graphs,  value  judgements  can  be  made  on 
goodness. 

The  graphs  of  the  six  time  integration  schemes  combined  with  the  linear  CFM  method  indicate  all 
of  these  techniques  to  be  diffusive  and  to  lag  behind  the  true  position  of  the  hill.  Figure  5.1.13  (c)  and 
(d)  contain  the  associated  Errors  for  the  linear  CFM.  First  note  that  the  order,  best  to  poorest,  is  the 
same  as  in  Figures  5.1.13  (a)  and  (b),  i.e,  standard,  modified,  and  smoothed.  From  the  Error  graph 
the  standard  Godunov  is  the  best  of  the  linear  CFM  methods  with  a  10%  Error  at  96  s.  If  the  Error 
graphs  for  linear  FEM  were  over-laid  on  those  for  linear  CFM,  all  of  the  respective  FEM  techniques 
would  have  less  error  than  the  associated  CFM  technique.  However,  the  difference  is  small  except  for 
the  standard  Lax-Wendroff.  Careful  observations  will  reveal  that  the  standard  finite  difference  tech¬ 
niques  (linear  CFM)  do  reasonably  well  compared  to  the  respective  linear  FEM  techniques. 

Whereas  the  cosine  hill  density  distribution  represents  a  "smooth"  function,  that  is  Equation  (2.2) 
and  its  first  derivative  are  continuous  functions,  the  square  hill  density  distribution  presents  a  more 
severe  test  for  the  solution  procedures. 
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The  square  hill  solutions  using  the  linear  FEM  are  presented  in  Figures  5.1.14  to  5.1.19  and  the 
associated  Error  graphs  in  Figures  5.1.26  (a)  and  (b).  After  observing  these,  it  is  noted  that  none  do  as 
well  as  they  did  for  the  cosine  hill.  The  smallest  Error  is  now  around  14%.  There  are  similarities 
between  the  two  propagated  hills.  The  two  standard  time  integration  schemes  still  track  the  hill  accu¬ 
rately.  The  standard  Lax-Wendroff  scheme  is  the  least  diffusive  (has  the  sharpest  front)  but  has  the 
largest  amplitude  oscillations  away  from  the  hill.  The  standard  Godunov  scheme  diffuses  the  hill  but  to 
a  large  extent  damps  out  the  oscillations.  From  observing  the  propagated  hill  graphs  and  the  Error 
graphs,  we  see  that  modified  Lax-Wendroff  and  Godunov  methods  give  results  comparable  to  the  stan¬ 
dard  methods.  The  smoothed  integrator  is  still  the  poorest  by  a  sizeable  margin. 

Turning  now  to  the  linear  CFM  method.  Figures  5.1.20  to  5.1.25  contain  results  of  the  advected 
square  hill  and  the  associated  Error  graphs  are  in  Figures  5.1.26(c)  and  (d).  First  note  that  the  pro¬ 
pagated  hill  does  not  look  any  worse  than  the  associated  linear  FEM  hill.  Next  note  from  the  Error 
graph  that  the  relative  positions  of  the  time  integration  schemes  are  the  same  as  previously  observed. 
The  smoothed  method  is  still  the  poorest  and  the  standard  and  the  modified  methods  give  nearly  the 
same  Error. 

Comparing  linear  FEM  and  linear  CFM  using  the  standard  and  modified  integrators  shows  that  at 
96  s  the  Errors  are  approximately  16%  (FEM)  and  22%  (CFM)  using  Lax-Wendroff  methods  and  15% 
(FEM  and  CFM)  using  Godunov  methods.  After  reviewing  the  associated  propagated  hill,  there  is  not 
much  difference  between  the  respective  FEM  and  CFM  methods.  For  example,  using  the  standard 
Lax-Wendroff  scheme,  the  propagated  hills  for  FEM  and  CFM  have  the  same  quality  of  tracking  the 
hill  and  maintaining  the  sharpness  of  the  front. 

5.2  Parabolic  Elements 

The  development  of  the  parabolic  finite  element  spatial  discretization  is  given  in  Section  3.2  and 
Appendix  B.  We  take  the  length  of  each  element  as  Z/  —  2  units  of  distance.  Thus  there  are  24  ele¬ 
ments  on  the  48  unit  long  track.  Recalling  that  each  parabolic  element  has  a  mid- point  node,  there  are 
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still  therefore  48  nodal  degrees  of  freedom.  By  using  24  parabolic  elements  a  direct  comparison  can  be 
made  with  48  linear  elements  on  the  basis  of  equal  "computational  effort  " 

The  advection  velocity  remains  at  V  =  1  unit  of  distance/s  and  the  time  step  remains  at  A i  =  0  2 
s  Ihe  Courant  number  based  on  element  length,  FA ill.',  is  0.10  for  the  parabolic  element  There  is 
a  question  as  whether  the  element  length  or  distance  between  node  points  is  the  appropriate  length 
scale  for  defining  the  Courant  number.  Based  on  the  latter  length  scale,  the  Courant  number  remains 
0.20 

The  graphs  for  the  smoothed  time  integration  schemes  are  not  included  for  the  parabolic  element 
f  rom  observations  discussed  in  Section  5.1,  it  was  shown  to  be  the  poorest  method  In  addition,  the 
smoothed  Lax-Wendroff  scheme  in  combination  with  the  parabolic  spatial  discretization  becomes 
unstable  in  a  short  time  period.  Error  values  greater  than  100"/.  were  found  in  less  than  24  s. 

Advection  of  the  cosine  hill  is  again  studied  first.  Figures  5.2.1  to  5.2.4  contain  the  parabolic 
FEM  solutions.  The  Error  graphs  are  in  Figures  5.2.9(a)  and  ( b ) .  Both  the  standard  and  modified 
Lax-Wendroff  time  integrations  do  an  excellent  job  of  tracking  the  hill  and  its  amplitude,  with  the 
modified  method  having  the  smallest  Error,  1.4%.  Both  Godunov  methods  track  the  hill  but  the  ampli¬ 
tude  is  reduced  by  numerical  diffusion.  The  Error  graphs  show  virtually  no  difference  between  the 
standard  and  modified  Godunov  methods. 

The  parabolic  CFM  solutions  for  the  advected  hill  are  shown  in  Figures  5  2  5  to  5  2  8.  The  asso¬ 
ciated  Error  graphs  appear  in  Figures  5.2.9(c)  and  (d)  For  the  two  Fax-Wendroff  methods,  parabolic 
CFM  does  remarkably  well,  giving  an  Error  of  2%  for  the  standard  and  4%  for  the  modified.  Also, 
both  Godunov  methods  give  results  nearly  identical  to  those  obtained  with  parabolic  FEM  Thus,  the 
condensing  of  the  mass  matrix  to  avoid  the  time  consuming  inversion  process  is  extremely  advan  a- 
geous  for  the  problem  of  the  cosine  hill 

Turning  now  to  the  advection  of  the  square  hill,  the  parabolic  FEM  results  appear  in  Figures 
5.2.10  to  5.2  13.  The  Error  graphs  are  in  Figures  5.2.18  (a)  and  (b).  As  with  linear  elements,  the 
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advected  solution  for  the  square  hill  is  not  as  accurate  as  that  for  the  cosine  hill.  All  the  parabolic  FEM 
schemes  track  the  hill  well.  The  standard  Lax-Wendroff  method,  which  gave  good  results  with  the 
linear  element  and  the  parabolic  element  for  a  cosine  hill  now  is  tending  to  loose  stability  at  96  s.  The 
modified  Lax-Wendroff  scheme  remains  stable  and  tracks  the  position  and  amplitude  of  the  hill  reason¬ 
ably  well.  The  Godunov  methods  retain  the  characteristics  of  diffusivity  and  positivity  that  have  been 
previously  observed. 

Figures  5.2.14  to  5.2.17  contain  results  of  the  parabolic  CFM  solutions  for  the  advected  square 
hill  The  Error  graphs  are  displayed  in  Figures  5.2.18  (c )  and  (d).  From  the  graphs  of  the  advected 
hills  there  does  not  seem  to  be  much  difference  in  the  tracking  abilities  of  either  CFM  or  FEM  com¬ 
bined  with  the  modified  Lax-Wendroff  scheme  or  either  Godunov  scheme.  It  is  worthwhile  to  note 
that,  for  the  standard  Lax-Wendroff  scheme,  the  CFM  is  stable  while  the  FEM  is  diverging  by  the  end 
of  96  s. 

After  examination  of  the  four  sets  of  Error  graphs  in  Figure  5.2.18.  it  is  seen  that  the  modified 
Lax-Wendroff  scheme  combined  with  parabolic  FEM  has  the  least  Error,  8%  at  96  s.  The  associated 
parabolic  CFM  has  an  Error  of  12%  at  96  s,not  a  significant  difference.  Upon  reviewing  the  advected 
square  hills  for  these  two  combinations,  no  significant  difference  is  found.  Both  combinations  track  the 
position  and  amplitude  of  the  hill  reasonably  well.  Both  produce  peak  amplitudes  that  are  approxi¬ 
mately  10%  greater  than  they  should  be  and  both  give  roughly  the  same  degree  of  sharpness  of  the  hill 
front. 

5.3  Comparison  of  Linear  and  Parabolic  Elements 

In  Table  5.3.1,  the  errors  generated  by  the  various  spatial  discretization/time  integration  schemes 
at  96  s  (two  cycles)  are  summarized.  The  "smoothed”  methods  are  not  included  in  the  table  since  their 
performances  were  markedly  inferior  to  the  other  methods  in  advecting  the  density  distribution. 

From  the  table,  we  see  that  the  numerically  diffusive,  standard  Godunov  time  integration  scheme 
overrides  whatever  spatial  discretization  is  used.  The  same  is  true  of  the  modified  Godunov  scheme 
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except  in  combination  with  linear  CFM.  For  the  linear  discretizations,  we  also  not  that  the  Godunov 
schemes  tend  to  outperform  the  Lax-Wendroff  schemes. 

Turning  our  attention  to  the  modified  Lax-Wendroff  time  integration  scheme,  we  not  diminishing 
errors  as  we  proceed  from  linear  CFM  to  linear  FEM  to  parabolic  CFM  to  parabolic  FEM.  This  latter 
combination  produces  the  least  overall  error  for  both  the  cosine  hill  and  square  hill  advection  problems. 

The  results  for  the  standard  Lax-Wendroff  scheme  are  mixed.  Linear  FEM  and  parabolic  CFM 
spatial  discretization  give  equivalent  errors  that  are  substantially  below  those  obtained  with  linear  CFM. 
However,  this  time  integration  in  combination  with  parabolic  FEM  performs  poorly  for  square  hill 
advection. 

5.4  Time  Step  Size  Study 

To  insure  that  the  results  of  the  previous  sections  were  not  biased  by  the  selection  of  A t  =  0.2  s 
and  also  to  examine  the  effects  of  time  step  size  on  the  stability  and  accuracy  of  the  advected  solution, 
a  study  of  time  step  variation  was  undertaken  for  some  of  the  better  performing  numerical  models. 

For  the  cosine  hill  density  distribution,  we  considered  the  standard  Lax-Wendroff,  linear  FEM 
combination  and  the  modified  Lax-Wendroff,  parabolic  FEM  combination.  The  Error  graphs  as  a  func¬ 
tion  of  time  step  size  are  given  in  Figure  5.4.1. 

For  the  standard  Lax-Wendroff,  linear  FEM  scheme,  we  see  that  the  error  decreases  as  A t 
decreases  from  0.5  s  to  0.1  s.  However,  there  is  no  significant  diffeicnce  between  the  errors  for 

At  =  0.2  s  and  At  =  0.1  s;  and,  in  fact,  an  asymptotic  error  of  around  1%  appears  to  have  been 

reached. 

In  contrast  to  the  decreasing  Error  with  decreasing  At,  the  modified  Lax-Wendroff,  parabolic  FEM 

scheme  appears  to  have  a  minimum  error  at  around  At  =*  0.25  s.  The  differences  in  error  for  time 

steps  ranging  from  At  -  0.1  s  to  At  =  0.33  s  are  not  meaningful,  however.  For  At  =  0.5  s,  this 
scheme  still  yields  an  acceptable  error  as  compared  to  the  previous  scheme  that  becomes  unstable.  The 
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advection  of  the  cosine  hill  with  Ar  =  0.5  s  is  shown  in  Figure  5.4.2.  The  computational  effort  here  is 
2/5th  of  that  required  for  Ar  =  0.2  s. 

Turning  now  to  the  square  hill,  the  effects  of  time  step  size  on  six  of  the  eight  combinations  of 
Lax-Wendroff  time  integration  and  spatial  discretization  were  examined.  The  Errors  for  linear  FEM  are 
shown  in  Figure  5.4.3,  and  for  linear  CFM  in  Figure  5.4.4.  We  see  from  these  figures  that  A r  =  0.2  s 
is  representative  of  the  general  behavior  of  each  scheme.  We  also  not  that,  as  with  the  cosine  hill,  the 
standard  Lax-Wendroff  scheme  gives  decreasing  error  with  decreasing  At,  while  the  modified  Lax- 
Wendroff  scheme  appears  to  have  a  minimum  error  at  an  intermediate  A/. 

The  Error  results  for  parabolic  FEM  are  given  in  Figure  5.4.5.  The  same  comments  as  made 
above  are  applicable  except  we  note  that  using  a  time  step  At  =  0.1  s  in  conjunction  with  standard 
Lax-Wendroff,  parabolic  FEM  provides  a  stable  solution. 

5.5  Time  Duration  Study 

As  a  final  check  on  the  performance  of  the  better  methods,  the  density  distributions  were  numeri¬ 
cally  advected  ten  times  (480  s)  around  the  racetrack.  Hence,  for  At  =  0.2  s,  there  were  2400  steps  in 
the  time  integration. 

Examining  the  results  for  the  cosine  hill  first.  Figures  5.5.1  and  5.5.2  show,  respectively,  the 
advected  solutions  for  the  standard  Lax-Wendroff,  linear  FEM  scheme  and  for  the  modified  Lax- 
Wendroff,  parabolic  FEM  scheme.  After  ten  cycles  around  the  track,  the  performance  of  the  latter 
scheme  is  seen  to  be  superior  to  the  performance  of  the  former  scheme.  The  Error  graph  for  this  prob¬ 
lem  is  shown  in  Figure  5.5.3.  A  4%  error  after  ten  cycles  around  the  track  is  exceptional  performance 
for  an  explicit  numerical  advector. 

For  the  square  hill,  the  modified  Lax-Wendroff  scheme  was  the  only  time  integrator  used.  The 
advected  solutions  for  this  time  integrator  in  combination  with  linear  FEM.  linear  CFM,  parabolic 
FEM,  and  parabolic  CFM  are  shown  in  Figures  5.5.4  to  5.5.7,  respectively.  After  ten  cycles,  it  is  seen 
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that  the  linear  FEM  and  linear  CFM  combinations  have  washed  out  most  of  the  information  contained 
in  the  original  square  hill.  The  parabolic  FEM  and  parabolic  C  FM  combinations,  meanwhile,  continue 
to  excellently  track  the  location  of  the  hill.  In  both  cases,  the  hill  has  lost  its  sharp  front  but  still  has  a 
distinct  front.  The  peak  amplitude  is  10%  greater  than  it  should  be.  Very  few  oscillations  occur  away 
from  the  hill.  For  both  cases,  the  Error  after  480  s  is  about  10%. 

Finally,  in  Figure  5.5  8,  the  advected  square  hill  for  the  modified  Lax-Wendroff,  parabolic  FEM 
scheme  with  A t  =  0.33  s  is  shown.  The  solution  is  again  quite  satisfactory.  The  Error  graph  comparing 
A /  =  0.2  s  and  A?  =  0.33  s  is  given  in  Figure  5.5.9. 

6.  CONCLUSIONS 

Numerical  studies  of  explicit  time  integration  techniques  for  the  finite  element  advection  equation 
have  been  conducted  The  results  have  shown  that,  overall,  the  "best”  numerical  solutions  are  obtained 
by  combining  the  modified  Lax-Wendroff  time  integration  with  parabolic  spatial  discretization  either  in 
its  lull  or  condensed  [ .V/ 1  matrix  form.  Both  the  standard  and  modified  Godunov  time  integrations 
have  been  shown  to  be  numerically  diffusive.  This  diffusivity  tends  to  override  whatever  spatial  discret¬ 
ization  is  used.  However,  the  positivity  property  possessed  by  the  Godunov  schemes  can  be  valuable 
for  certain  applications. 

Many  numerical  questions  remain  for  future  studies  to  resolve.  These  include  the  abilities  of  the 
belter  schemes  to  integrate  coupled  sets  of  nonlinear  hyperbolic  equations  and  to  handle  multidimen¬ 
sional  problems  This  latter  issue  can  grow  exponentially  since  we  have  shown  the  better  schemes  to 
involve  nonlinear  shape  functions.  Many  analytical  questions  regarding  the  stability  and  phase  proper¬ 
ties  of  the  various  schemes  also  need  to  be  answered 
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Appendix  A 

LINEAR  ELEMENT  MATRICES 

Let  the  x  domain  be  subdivided  into  N  linear  finite  elements  each  of  length  U.  The  nodes  are 
numbered  consecutively  from  1  to  N  +  1  as  x  increases  from  0  to  L.  The  element  mass  [A/']  and 
advection  [AT*]  matrices  are  given  by  Equations  (3.9)  and  (3.10). 

On  recalling  that  nodes  1  and  N  +  1  represent  the  same  point  (racetrack  boundary  condition),  the 
assembled  mass  [M]  and  advection  (A]  matrices  for  the  A  degree  of  freedom  system  are  found  as 
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Appendix  B 

PARABOLIC  ELEMENT  MATRICES 

Lei  a  domain  be  subdivided  into  S  parabolic  finite  elements  each  of  length  L‘ .  The  nodes  are 
numbered  consecutively  from  1  to  2A  +  1  as  v  increase  from  0  to  L.  The  even  numbered  nodes 
correspond  to  element  mid-point  nodes.  The  element  mass  |.V/'J  and  advection  [A"']  matrices  are  given 
by  equations  (3.13)  and  (3.14). 

On  recalling  that  nodes  1  and  2 N  +  1  represent  the  same  point  (racetrack  boundary  condition), 
the  assembled  mass  [M]  and  advection  [A]  matrices  for  the  IS  degree  of  freedom  system  are  found  as 
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Fig.  5.1.25  —  Advectcd  square  hill:  modified  Godunov,  linear  CFM 
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Fig.  5.1.26  —  Average  absolute  errors  for  square  hill  advection  (a)  linear  FEM,  Lax-Wcndrofl  time  integrators  tbl 
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Fig  5.1.26  (Continued)  —  Average  absolute  errors  for  square  hill  advection.  (a)  linear  F1M.  Lax-W'cndrolT  time 
integrators  (b)  linear  FEM.  Godunov  time  integrators  (c)  linear  C'FM,  Lax-WcndrofT  time  integrators  (d)  linear  C'FM. 
Godunov  time  integrators 
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Fig.  5.2.3  —  Advected  cosine  hill:  standard  Godunov,  parabolic  FEM 
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Fig.  5.2.6  —  Advected  cosine  hill:  modified  Lax-Wendroff, 
parabolic  CFM 
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Fig  5.2  12  —  Advecled  square  hill:  standard  Godunov, 
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Fig.  5.2.17  —  Advected  square  hill:  modified  Godunov, 
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Fig  5.2.18  —  Average  absolute  errors  for  square  hill  advection  (a)  parabolic  FEM,  Lax-Wendroff  time  integrators  (b) 
parabolic  FEM,  Godunov  time  integrators  (c)  parabolic  CFM,  Lax-Wendroff  time  integrators  (d)  parabolic  CFM, 
Godunov  time  integrators 
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Fig.  5.4.3  —  Average  absolute  errors  for  square  hill  advection  versus  time  step  size  using  linear  FEM 
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